!! Copyright (C) 2009,2010,2011,2012  Marco Restelli
!!
!! This file is part of:
!!   LDGH -- Local Hybridizable Discontinuous Galerkin toolkit
!!
!! LDGH is free software: you can redistribute it and/or modify it
!! under the terms of the GNU General Public License as published by
!! the Free Software Foundation, either version 3 of the License, or
!! (at your option) any later version.
!!
!! LDGH is distributed in the hope that it will be useful, but WITHOUT
!! ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
!! or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
!! License for more details.
!!
!! You should have received a copy of the GNU General Public License
!! along with LDGH. If not, see <http://www.gnu.org/licenses/>.
!!
!! author: Marco Restelli                   <marco.restelli@gmail.com>


!>\brief
!!
!! Collection of multistep time integrators
!!
!! \n
!!
!! The methods implemented in this module are aimed at solving the ODE
!! problem
!! \f{displaymath}{
!!  u' = f(t,u).
!! \f}
!!
!! A generic \f$p+1\f$-step method has the form
!! \f{displaymath}{
!!  u_{n+1} = \sum_{j=0}^p a_ju_{n-j} + \Delta t
!!  \sum_{j=-1}^pb_jf_{n-j}.
!! \f}
!! The method is of order \f$q\f$ if it satisfies the following
!! compatibility conditions:
!! \f{displaymath}{
!!  \sum_{j=0}^p (-j)^ia_j + i \sum_{j=-1}^p(-j)^{i-1}b_j = 1, \qquad
!!  i=0,\ldots,q.
!! \f}
!!
!! The initialization of a multistep method must be coherent with the
!! order of the method. In particular, to initialize a method of order
!! \f$q\f$, we have to ensure an error \f$\mathcal{O}( \Delta t^q )\f$
!! on each starting value. Typically, this implies the use of a method
!! of order \f$q-1\f$.
!!
!! \warning Due to the initialization error, the numerical solution of
!! the test problem
!! \f{displaymath}{
!!   u' = t^{q-1}
!! \f}
!! will not be exact. However, the error should remain constant during
!! the integration.
!! 
!! \note The initialization subroutines of the multistep methods also
!! perform the first \f$p\f$ steps; this must be taken into account to
!! keep track of the sistem time.
!!
!! \section BDF BDF schemes
!!
!! BDF schemes are characterized by the fact that \f$b_{-1}\neq0\f$
!! while \f$b_{j}=0\f$ for \f$j\geq0\f$. This implies that the right
!! hand side \f$f\f$ of the problem is never computed explicitly, but
!! only as part of the solution of the implicit problem
!! \f{displaymath}{
!!   u - \Delta t\, b_{-1}\, f(t,u) = b.
!! \f}
!! In practical terms, this means that the member function \c rhs of
!! the class \c c_ode is not used by BDF methods. With this respect,
!! however, there are two issues that deserve attention:
!! <ul>
!!  <li> the computation of the linearization state \f$x_l\f$ in
!!  principle could make explicit use of the right hand side of the
!!  equation \f$f\f$; in the present module, all the BDF schemes
!!  compute \f$x_l\f$ with an extrapolation of the last time steps,
!!  without making explicit use of the right hand side;
!!  <li> the initialization steps in principle could make explicit use
!!  of the right hand side of the equation; again, all the BDF schemes
!!  implemented in this module use BDF type schemes of lower order as
!!  initialization schemes.
!! </ul>
!!
!! \section constrained Constrained problems
!!
!! Some of the solvers defined in the present module can also be used
!! to integrate problems of the form
!! \f{displaymath}{
!!  \begin{array}{rcl}
!!   u' & = & f(t,u,p) \\
!!   g(t,u) & = & 0
!!  \end{array}
!! \f}
!! fore some constraint function \f$g\f$ and associated Lagrange
!! multiplier \f$p\f$. Because of the constraint, an
!! implicit method will be required.
!!
!! Using BDF for constrained problems is straightforward, essentially
!! because such solvers do not make any use of the member function \c
!! rhs of the ODE class \c c_ode. In fact, the only required
!! modification is implementing the \c solve member function so that
!! it solves
!! \f{displaymath}{
!!  \begin{array}{rcl}
!!   u - \sigma f(t,u,p) & = & b_u \\
!!   g(t,u) & = & 0
!!  \end{array}
!! \f}
!! for arbitrary \f$\sigma\f$ and \f$b_u\f$, where \f$b_u\f$ denotes
!! the components of \f$b\f$ corresponding to the prognostic equations
!! for \f$u\f$ (the remaining components of the right hand side
!! \f$b\f$ can be ignored in \c solve).
!<----------------------------------------------------------------------
module mod_multistep

!-----------------------------------------------------------------------

 use mod_messages, only: &
   mod_messages_initialized, &
   error,   &
   warning, &
   info

 use mod_kinds, only: &
   mod_kinds_initialized, &
   wp

 use mod_time_integrators_base, only: &
   mod_time_integrators_base_initialized, &
   c_ode

!-----------------------------------------------------------------------
 
 implicit none

!-----------------------------------------------------------------------

! Module interface

 public :: &
   mod_multistep_constructor, &
   mod_multistep_destructor,  &
   mod_multistep_initialized, &
   thetam_init, thetam_step, thetam_clean, &
   bdf1_init,   bdf1_step,   bdf1_clean,   &
   bdf2_init,   bdf2_step,   bdf2_clean,   &
   bdf3_init,   bdf3_step,   bdf3_clean,   &
   bdf2ex_init, bdf2ex_step, bdf2ex_clean

 private

!-----------------------------------------------------------------------

! Module variables

 ! public members
 logical, protected ::               &
   mod_multistep_initialized = .false.

 ! private members

 !> time step
 real(wp) :: dt
 !> theta parameter
 real(wp) :: th

 character(len=*), parameter :: &
   this_mod_name = 'mod_multistep'

!-----------------------------------------------------------------------

contains

!-----------------------------------------------------------------------

 subroutine mod_multistep_constructor()
  character(len=*), parameter :: &
    this_sub_name = 'constructor'

   !Consistency checks ---------------------------
   if( (mod_messages_initialized             .eqv..false.) .or. &
       (mod_kinds_initialized                .eqv..false.) .or. &
       (mod_time_integrators_base_initialized.eqv..false.) ) then
     call error(this_sub_name,this_mod_name, &
                'Not all the required modules are initialized.')
   endif
   if(mod_multistep_initialized.eqv..true.) then
     call warning(this_sub_name,this_mod_name, &
                  'Module is already initialized.')
   endif
   !----------------------------------------------

   mod_multistep_initialized = .true.
 end subroutine mod_multistep_constructor

!-----------------------------------------------------------------------
 
 subroutine mod_multistep_destructor()
  character(len=*), parameter :: &
    this_sub_name = 'destructor'
   
   !Consistency checks ---------------------------
   if(mod_multistep_initialized.eqv..false.) then
     call error(this_sub_name,this_mod_name, &
                'This module is not initialized.')
   endif
   !----------------------------------------------

   mod_multistep_initialized = .false.
 end subroutine mod_multistep_destructor

!-----------------------------------------------------------------------
 
 !> \f$\theta\f$-method initialization
 subroutine thetam_init(dt_step,theta)
  real(wp), intent(in) :: dt_step, theta
 
  character(len=*), parameter :: &
    this_sub_name = 'thetam_init'

   dt = dt_step
   th = theta
 
 end subroutine thetam_init
 
!-----------------------------------------------------------------------
 
 !> \f$\theta\f$-method step
 !!
 !! The method is
 !! \f{displaymath}{
 !!  u_{n+1} - \Delta t \, \theta f(u_{n+1}) = u_n + \Delta t \,
 !!  (1-\theta) f(u_{n}).
 !! \f}
 !!
 !! In practice, we allow for an approximation (typically, a
 !! linearization) of the implicit term, so that the method becomes
 !! \f{displaymath}{
 !!  u_{n+1} - \Delta t \, \theta \tilde{f}_{x_l}(u_{n+1}) = u_n +
 !!  \Delta t \, (1-\theta) f(u_{n}).
 !! \f}
 !! Concerning the state used for the linearization \f$x_l\f$, we
 !! notice that, in order to obtain a second order method when
 !! \f$\theta=1/2\f$, we need at least a second order approximation of
 !! \f$x_{n+1}\f$. In fact, making the substitution
 !! \f{displaymath}{
 !!  f(t_{n+1},u_{n+1}) \to \tilde{f} + \Delta f,
 !! \f}
 !! we obtain the local trucation error
 !! \f{displaymath}{
 !!  \begin{array}{rcl}
 !!   \displaystyle
 !!   \Delta t \, \tau(\Delta t) & = &
 !!   \displaystyle
 !!   y_{n+1}-y_n-\frac{\Delta t}{2}\left[ f(t_{n+1},y_{n+1}) - \Delta
 !!   f +f(t_{n},y_{n}) \right] \\[3mm]
 !!   & = &
 !!   \displaystyle
 !!   \frac{\Delta t}{2}\Delta f + \mathcal{O}(\Delta t^3).
 !!  \end{array}
 !! \f}
 !! For \f$f\f$ smooth, a second order approximation \f$x_l\f$ of
 !! \f$y_{n+1}\f$ guarantees \f$\Delta f=\mathcal{O}(\Delta t^2)\f$,
 !! preserving the order of the scheme.
 !<
 subroutine thetam_step(unew,t,unow,tmp)
  real(wp),     intent(in)    :: t
  class(c_ode), intent(in)    :: unow
  class(c_ode), intent(inout) :: unew
  !> \todo this argument should become a module variable
  class(c_ode), intent(inout) :: tmp(:)
 
  character(len=*), parameter :: &
    this_sub_name = 'thetam_step'

   ! compute the tendency
   call tmp(1)%rhs(t,unow)
   ! compute the rhs
   call tmp(2)%op( unow , (1.0_wp-th)*dt , tmp(1) , unew )
   ! compute the linearization state by explicit Euler
   call tmp(3)%op( unow ,      dt        , tmp(1) , unew )
   ! solve the implicit problem
   call unew%solve( t+dt , th*dt , tmp(2) , tmp(3) )
 
 end subroutine thetam_step
 
!-----------------------------------------------------------------------
 
 !> \f$\theta\f$-method clean-up
 pure subroutine thetam_clean()
  character(len=*), parameter :: &
    this_sub_name = 'thetam_clean'

  ! nothing to do
 
 end subroutine thetam_clean
 
!-----------------------------------------------------------------------

 !> BDF1 initialization
 !!
 !! Nothing to do, indeed.
 subroutine bdf1_init(dt_step)
  real(wp),     intent(in)    :: dt_step
 
  character(len=*), parameter :: &
    this_sub_name = 'bdf1_init'

   dt = dt_step

 end subroutine bdf1_init
 
!-----------------------------------------------------------------------
 
 !> BDF1 step
 !!
 !! This is inpractice the implicit Euler method
 !! \f{displaymath}{
 !!  u_{n+1} - \Delta t\, f(u_{n+1}) = u_n.
 !! \f}
 !! This method differs from the \f$\theta\f$-method with
 !! \f$\theta=1\f$ because we use here \f$x_l=u_n\f$, i.e. there is
 !! extrapolation for the linearization state.
 pure subroutine bdf1_step(unew,t,unow)
  real(wp),     intent(in)    :: t
  class(c_ode), intent(in)    :: unow
  class(c_ode), intent(inout) :: unew
 
  character(len=*), parameter :: &
    this_sub_name = 'bdf1_step'

   call unew%solve( t+dt , dt , unow , unow )

 end subroutine bdf1_step
 
!-----------------------------------------------------------------------
 
 !> BDF1 clean-up
 pure subroutine bdf1_clean()
  character(len=*), parameter :: &
    this_sub_name = 'bdf1_clean'

  ! nothing to do
 
 end subroutine bdf1_clean
 
!-----------------------------------------------------------------------

 !> BDF2 initialization
 !!
 !! This subroutine performs the first time step using the implicit
 !! Euler method and copies \c unow in \c uold.
 subroutine bdf2_init(dt_step,t,unow,uold)
  real(wp),     intent(in)    :: dt_step, t
  class(c_ode), intent(inout) :: unow
  !> \todo this arguments should become allocatable module variables
  class(c_ode), intent(inout) :: uold
 
  character(len=*), parameter :: &
    this_sub_name = 'bdf2_init'

   dt = dt_step

   ! set uold
   call uold%assign(unow)

   ! compute the first time step (implicit Euler)
   call unow%solve( t+dt , dt , uold , uold )

 end subroutine bdf2_init
 
!-----------------------------------------------------------------------
 
 !> BDF2 step
 !!
 !! The method is
 !! \f{displaymath}{
 !!  u_{n+1} - \Delta t \, b_{-1}f(u_{n+1}) = a_0u_n + a_1u_{n-1}
 !! \f}
 !! with \f$a_0=4/3\f$, \f$a_1=-1/3\f$ and \f$b_{-1}=2/3\f$.
 !!
 !! In practice, we allow for an approximation (typically, a
 !! linearization) of the implicit term, so that the method becomes
 !! \f{displaymath}{
 !!  u_{n+1} - \Delta t \, b_{-1}\tilde{f}_{x_l}(u_{n+1}) = a_0u_n +
 !!  a_1u_{n-1}.
 !! \f}
 !! To ensure second order accuracy, the linearization state \f$x_l\f$
 !! is extrapolated using the last two time steps:
 !! \f{displaymath}{
 !!  x_l = 2u_n - u_{n-1}.
 !! \f}
 !<
 pure subroutine bdf2_step(unew,t,unow,tmp,uold)
  real(wp),     intent(in)    :: t
  class(c_ode), intent(in)    :: unow
  class(c_ode), intent(inout) :: unew
  !> \todo this argument should become a module variable
  class(c_ode), intent(inout) :: tmp(:), uold
 
  real(wp), parameter :: &
    a0 =  4.0_wp/3.0_wp, &
    a1 = -1.0_wp/3.0_wp, &
    bm =  2.0_wp/3.0_wp
  character(len=*), parameter :: &
    this_sub_name = 'bdf2_step'

   ! compute the rhs of the implicit problem and the extrapolation 
   call tmp(1)%stimes(a1,uold)
   call unew%stimes(-1.0_wp,uold)        ! uold not needed any more
   call tmp(2)%op(tmp(1),a0,unow,uold)   ! rhs
   call tmp(1)%op(unew,2.0_wp,unow,uold) ! extrapolation
   ! solve the implicit problem
   call unew%solve( t+dt , bm*dt , tmp(2) , tmp(1) )

   ! update the multistep history
   call uold%assign(unow)

 end subroutine bdf2_step
 
!-----------------------------------------------------------------------
 
 !> BDF2 clean-up
 pure subroutine bdf2_clean()
  character(len=*), parameter :: &
    this_sub_name = 'bdf2_clean'

  ! nothing to do
 
 end subroutine bdf2_clean
 
!-----------------------------------------------------------------------

 !> BDF3 initialization
 !!
 !! This subroutine performs the first two steps using the BDF2 method
 !! and copies \c unow in \c uold.
 subroutine bdf3_init(dt_step,t,unow,tmp,uold)
  real(wp),     intent(in)    :: dt_step, t
  class(c_ode), intent(inout) :: unow
  !> \todo this argument should become an allocatable module variable
  class(c_ode), intent(inout) :: tmp(:), uold(:)
 
  real(wp), parameter :: &
    a0 =  4.0_wp/3.0_wp, &
    a1 = -1.0_wp/3.0_wp, &
    bm =  2.0_wp/3.0_wp
  character(len=*), parameter :: &
    this_sub_name = 'bdf3_init'

   dt = dt_step

   ! set uold(2)
   call uold(2)%assign(unow)

   ! leap-frog type step for the first step: set uold(1)
   !
   ! We solve the problem
   !     (u1-u0)/dt = f((u1+u0)/2)
   ! defining   z = (u1+u0)/2. To obtain a good linearization state we
   ! make one iteration.
   call unow%solve( t+dt/2.0_wp , dt/2.0_wp , uold(2) , uold(2) )
   call uold(1)%solve( t+dt/2.0_wp , dt/2.0_wp , uold(2) , unow )
   call tmp(1)%stimes(2.0_wp,uold(1)) ! compute u1 = 2*z - u0
   call uold(1)%op(tmp(1),-1.0_wp,uold(2),tmp(3))

   ! BDF2 step: set unow
   call unow%stimes(   a1  ,uold(2))
   call tmp(1)%op(unow,  a0  ,uold(1),tmp(3)) ! rhs
   call unow%stimes(-1.0_wp,uold(2))
   call tmp(2)%op(unow,2.0_wp,uold(1),tmp(3)) ! extrapolation
   call unow%solve( t+2.0_wp*dt , bm*dt , tmp(1) , tmp(2) )

   ! Notice that this initialization requires three tmp, however for
   ! the time stepping 2 tmp are enough.

 end subroutine bdf3_init
 
!-----------------------------------------------------------------------
 
 !> BDF3 step
 !!
 !! The method is
 !! \f{displaymath}{
 !!  u_{n+1} - \Delta t \, b_{-1}f(u_{n+1}) = a_0u_n + a_1u_{n-1} +
 !!  a_2u_{n-2}
 !! \f}
 !! with \f$a_0=18/11\f$, \f$a_1=-9/11\f$, \f$a_2=2/11\f$, and
 !! \f$b_{-1}=6/11\f$.
 !!
 !! To ensure third order accuracy, the linearization state \f$x_l\f$
 !! is extrapolated using the last three time steps:
 !! \f{displaymath}{
 !!  x_l = 3u_n - 3u_{n-1} + u_{n-2}.
 !! \f}
 !<
 pure subroutine bdf3_step(unew,t,unow,tmp,uold)
  real(wp),     intent(in)    :: t
  class(c_ode), intent(in)    :: unow
  class(c_ode), intent(inout) :: unew
  !> \todo this argument should become a module variable
  class(c_ode), intent(inout) :: tmp(:), uold(:)
 
  real(wp), parameter :: &
    a0 =  18.0_wp/11.0_wp, &
    a1 = - 9.0_wp/11.0_wp, &
    a2 =   2.0_wp/11.0_wp, &
    bm =   6.0_wp/11.0_wp
  character(len=*), parameter :: &
    this_sub_name = 'bdf3_step'

   ! compute the rhs of the implicit problem and the extrapolation
   call tmp(1)%stimes(a2,uold(2))
   call unew%op(uold(2),-3.0_wp,uold(1),tmp(2)) ! done with uold(2)
   call tmp(2)%op(tmp(1),a1,uold(1),uold(2))
   call tmp(1)%op(tmp(2),a0,unow,uold(2))
   call tmp(2)%op(unew,3.0_wp,unow,uold(2))
   ! solve the implicit problem
   call unew%solve( t+dt , bm*dt , tmp(1) , tmp(2) )

   ! update the multistep history
   call uold(2)%assign(uold(1))
   call uold(1)%assign(unow)

 end subroutine bdf3_step
 
!-----------------------------------------------------------------------
 
 !> BDF2 clean-up
 pure subroutine bdf3_clean()
  character(len=*), parameter :: &
    this_sub_name = 'bdf3_clean'

  ! nothing to do
 
 end subroutine bdf3_clean
 
!-----------------------------------------------------------------------

 !> BDF2ex initialization
 subroutine bdf2ex_init(dt_step,t,unow,tmp,uold)
  real(wp),     intent(in)    :: dt_step, t
  class(c_ode), intent(inout) :: unow
  class(c_ode), intent(inout) :: tmp(:), uold(:)
 
  character(len=*), parameter :: &
    this_sub_name = 'bdf2ex_init'

   dt = dt_step

   call uold(1)%assign(unow)

   call tmp(1)%rhs(t,unow)
   call unow%op( uold(1) , dt , tmp(1) , uold(2) )

   call uold(2)%rhs(t,uold(1))
 
 end subroutine bdf2ex_init
 
!-----------------------------------------------------------------------
 
 !> BDF2ex step
 !!
 !! The method is similar to BDF2, except that we use an extrapolation
 !! for the right hand side at time \f$n+1\f$ so that the resulting
 !! method is explicit. In practice, we have
 !! \f{displaymath}{
 !!  u_{n+1} = a_0u_n + a_1u_{n-1} + \Delta t \left[ \, b_0 f(u_n) +
 !!  b_1 f(u_{n-1}) \right],
 !! \f}
 !! with \f$a_0=4/3\f$, \f$a_1=-1/3\f$, \f$b_{0}=4/3\f$ and
 !! \f$b_{1}=-2/3\f$.
 !<
 subroutine bdf2ex_step(unew,t,unow,tmp,uold)
  real(wp),     intent(in)    :: t
  class(c_ode), intent(in)    :: unow
  class(c_ode), intent(inout) :: unew
  !> \todo this argument should become a module variable
  class(c_ode), intent(inout) :: tmp(:), uold(:)
 
  real(wp), parameter :: &
    a0 =  4.0_wp/3.0_wp, &
    a1 = -1.0_wp/3.0_wp, &
    b0 =  4.0_wp/3.0_wp, &
    b1 = -2.0_wp/3.0_wp
  character(len=*), parameter :: &
    this_sub_name = 'bdf2ex_step'

   ! compute the "a" terms
   call tmp(1)%stimes(a1,uold(1))
   call tmp(2)%op(tmp(1),a0,unow,uold(1))
   ! compute the "b" terms
   call uold(1)%stimes(b1,uold(2))
   call uold(2)%rhs(t,unow)
   call tmp(1)%op(uold(1),b0,uold(2),unew)
   ! sum the terms
   call unew%op(tmp(2),dt,tmp(1),uold(1))

   ! update the multistep history
   call uold(1)%assign(unow)
   ! uold(2) already has the correct value

 end subroutine bdf2ex_step
 
!-----------------------------------------------------------------------
 
 !> BDF2ex clean-up
 pure subroutine bdf2ex_clean()
  character(len=*), parameter :: &
    this_sub_name = 'bdf2ex_clean'

  ! nothing to do
 
 end subroutine bdf2ex_clean
 
!-----------------------------------------------------------------------

end module mod_multistep
